;; Non-robust matrix functions
;; Copyright (C) 2021  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.linalg.matrix)

;; A matrix is defined as a map: {:ncols ncols, :data data}, where
;; ncols is the number of columns in the matrix;
;; data is a sequence of matrix data, in row-major order.
;;
;; For example, the matrix:
;; 0 1 2 3
;; 4 5 6 7
;; can be represented by:
;; {:ncols 4 :data [0 1 2 3 4 5 6 7]}
;; or even:
;; {:ncols 4 :data (range 8)}
;;
;; If the number of elements in data is not a multiple of ncols then the
;; matrix is undefined and the unexpected probably will happen.

(defn- index
  "Convert row `i` column `j` to an index to the matrix data with `ncols` columns."
  [i j ncols]
  (+ (* i ncols) j))

(defn- end-row
  "Index of the end of row `i`, exclusive."
  [i ncols]
  (* (inc i) ncols))

(defn column
  "Return column `j` beginning at row `i`."
  ([i j {:keys [ncols data]}]
   (column i j ncols data))
  ([i j ncols data]
   (->> data
        (drop (index i j ncols))
        (take-nth ncols))))

(defn columns
  "Return a sequence of columns, each column being a sequence of elements."
  ([{:keys [ncols data]}]
   (columns ncols data))
  ([ncols data]
   (letfn [(f [i ncols data]
             (if (< i ncols)
               (lazy-seq (cons (take-nth ncols data) (f (inc i) ncols (rest data))))
               []))]
     (f 0 ncols (seq data)))))

(defn columns->matrix
  "Create a matrix from a sequence of columns.
  All columns must have the same number of elements otherwise the
  result is undefined."
  [columns]
  {:ncols (count columns) :data (apply interleave columns)})

(defn row
  "Return row `i` beginning at column `j`."
  ([i j {:keys [ncols data]}]
   (row i j ncols data))
  ([i j ncols data]
   (subvec data (index i j ncols) (end-row i ncols))))

(defn rows
  "Return a sequence of rows, each row being a sequence."
  ([{:keys [ncols data]}]
   (rows ncols data))
  ([ncols data]
   (partition ncols data)))

(defn rows->matrix
  "Create a matrix from a sequence of rows.
  All rows must have the same number of elements otherwise the result
  is undefined."
  [rows]
  {:ncols (count (first rows)) :data (apply concat rows)})
